Analisis Descriptivo de Series de Tiempo en R
Gráfico de una Serie de tiempo
Vamos a analizar de forma descriptiva algunas serie de tiempo. Empezaremos por la serie de desempleo de los Estados Unidos que viene en el paquete TSstudio.
library(TSstudio)
data(USUnRate)
ts_info(USUnRate)## The USUnRate series is a ts object with 1 variable and 864 observations
## Frequency: 12
## Start time: 1948 1
## End time: 2019 12class(USUnRate)## [1] "ts"plot(USUnRate,main = "US Monthly Unemployment Rate",ylab="Unemployment Rate (%)")
La serie es medida mensual, es decir, presenta una frecuencia de 12. Qué
características podemos observar? * Tendencia? * Estacionalidad? *
Cíclos? * Varianza marginal no constante? Vamos a seleccionar un periodo
de tiempo mas corto
unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
title = "US Monthly Unemployment Rate",
Ytitle = "Unemployment Rate (%)",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)Note que aquí podemos ver varias características: Estacionalidad(no es tan evidente), tres periodos de ciclos, el primero de 1990 a 2000,el segundo de 2000 a 2007, y el tercero de 2007 a 2019. No parece tener una heterocedasticidad marginal.
Veamos ahora la tasa de desempleo de Colombia. Hay que hacerle un ajuste a la base de datos porque está en orden descendente en el tiempo.
Desempleo y Empleo Colombia
library(readxl)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()DesempleoyEmpleo <- read_excel("DesempleoyEmpleo.xlsx", range="A9:C249")
str(DesempleoyEmpleo)## tibble [240 × 3] (S3: tbl_df/tbl/data.frame)
## $ Fecha : chr [1:240] "2020-12" "2020-11" "2020-10" "2020-09" ...
## $ Tasadeempleo : num [1:240] 53.4 53.2 53.2 50.6 49.3 ...
## $ Tasadedesempleo: num [1:240] 13.4 13.3 14.7 15.8 16.8 ...DesempleoyEmpleo_1=DesempleoyEmpleo %>% map_df(rev)
tail(DesempleoyEmpleo)## # A tibble: 6 × 3
## Fecha Tasadeempleo Tasadedesempleo
## <chr> <dbl> <dbl>
## 1 2001-06 51.8 15.2
## 2 2001-05 51.2 14.2
## 3 2001-04 51.5 14.6
## 4 2001-03 53.0 15.7
## 5 2001-02 52.7 17.3
## 6 2001-01 53.0 16.7head(DesempleoyEmpleo_1)## # A tibble: 6 × 3
## Fecha Tasadeempleo Tasadedesempleo
## <chr> <dbl> <dbl>
## 1 2001-01 53.0 16.7
## 2 2001-02 52.7 17.3
## 3 2001-03 53.0 15.7
## 4 2001-04 51.5 14.6
## 5 2001-05 51.2 14.2
## 6 2001-06 51.8 15.2library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(xts)##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, lastFechas=as.yearmon(DesempleoyEmpleo_1$Fecha)
Desempleo_Col_xts=xts(x = DesempleoyEmpleo_1$Tasadedesempleo,frequency = 12,order.by = Fechas)
ts_info(Desempleo_Col_xts)## The Desempleo_Col_xts series is a xts object with 1 variable and 240 observations
## Frequency: monthly
## Start time: Jan 2001
## End time: Dec 2020plot(Desempleo_Col_xts)
ts_plot(Desempleo_Col_xts,
title = "Tasa de Desemplo Mensual Colombia",
Ytitle = "Tasa de Desempleo(%)",
Xtitle = "Año",
Xgrid = TRUE,
Ygrid = TRUE)Qué características presenta esta serie?
Otros objetos
Note que también se puede crear un objeto tsibble y usar el paquete feast, fable y fabletools(tsibble) y timetk(tibble), lo cuales funcionan con el paquete tidyverse.
require(feasts)## Loading required package: feasts## Loading required package: fabletoolsrequire(fable)## Loading required package: fablerequire(timetk)## Loading required package: timetkrequire(tsibble)## Loading required package: tsibble##
## Attaching package: 'tsibble'## The following object is masked from 'package:zoo':
##
## index## The following objects are masked from 'package:base':
##
## intersect, setdiff, union###Creación objeto tssible a partir de un objeto tibble
df_desempleo=data.frame(Desempleo=DesempleoyEmpleo_1$Tasadedesempleo,Fecha=Fechas)
tbl_desempleo=tibble(df_desempleo)
tbl_desempleo_format_fecha=tbl_desempleo
tbl_desempleo_format_fecha$Fecha=yearmonth(tbl_desempleo_format_fecha$Fecha)
###El tipo de fechas debe ser alguno que reconozca tsibble
tsbl_desempleo=as_tsibble(tbl_desempleo_format_fecha,index=Fecha) ####La fecha en tsibble es importante
##Gráfica de tsibble
autoplot(tsbl_desempleo,Desempleo)+labs(tittle="Serie de Desempleo Colombia Mensual",y="Tasa de Deszempleo")
###Gráfica timetk
tbl_desempleo%>%plot_time_series(.value=Desempleo,.date_var=Fecha)Análisis de Tendencias
Vamos a ver la forma de estimar la tendencia y/o eliminarla.
library(astsa)
library(TSstudio)
data(chicken)
ts_info(chicken)## The chicken series is a ts object with 1 variable and 180 observations
## Frequency: 12
## Start time: 2001 8
## End time: 2016 7plot(chicken,main="Precio Mensual de la Libra de Pollo en Estados Unidos", ylab="Precio en Centavos de Dólar")
#ts_plot(chicken)Al parecer la serie de precios mensuales del pollo presenta una tendencia creciente al parecer lineal, es decir \[y_t=\mu_t+a_t\] o mas específicamente
\[y_t=\beta_0+\beta_1 t +a_t\]
summary(fit <- lm(chicken~time(chicken), na.action=NULL))##
## Call:
## lm(formula = chicken ~ time(chicken), na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7411 -3.4730 0.8251 2.7738 11.5804
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.131e+03 1.624e+02 -43.91 <2e-16 ***
## time(chicken) 3.592e+00 8.084e-02 44.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.696 on 178 degrees of freedom
## Multiple R-squared: 0.9173, Adjusted R-squared: 0.9168
## F-statistic: 1974 on 1 and 178 DF, p-value: < 2.2e-16plot(chicken, ylab="centavos por libra")
abline(fit,col = "red") # Se añade la recta ajusta
###Eliminamos la tendencia con la predicción la recta
ElimiTendchick=chicken-predict(fit)
plot(ElimiTendchick,main="Serie Chicken Sin tendencia")
library(tidyverse)
library(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:tsibble':
##
## interval## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(timetk)
library(tsibble)
# Se configura para gráficas plotly(# FALSE retorna ggplots y no plotly)
interactive <- FALSE
indice_chicken=as.Date(as.yearmon(tk_index(chicken)))
## Otra forma de extraer el indice estimetk::tk_index(chicken)
df_chicken=data.frame(Fecha=indice_chicken,Pollo=as.matrix(chicken))
str(df_chicken)## 'data.frame': 180 obs. of 2 variables:
## $ Fecha: Date, format: "2001-08-01" "2001-09-01" ...
## $ Pollo: num 65.6 66.5 65.7 64.3 63.2 ...tibble_chicken=tibble(df_chicken)
duplicates(tibble_chicken, key = NULL, index=Fecha) ##Mirar si hay registros duplicados## # A tibble: 0 × 2
## # … with 2 variables: Fecha <date>, Pollo <dbl>
## # ℹ Use `colnames()` to see all variable namestibble_chicken_fechas_correct=tibble_chicken
tibble_chicken_fechas_correct$Fecha=yearmonth(tibble_chicken_fechas_correct$Fecha)
print(duplicates(tibble_chicken_fechas_correct, key = NULL, index=Fecha))## # A tibble: 0 × 2
## # … with 2 variables: Fecha <mth>, Pollo <dbl>
## # ℹ Use `colnames()` to see all variable namestsibble_chicken=tsibble(tibble_chicken_fechas_correct,index=Fecha)
tsibble_chicken## # A tsibble: 180 x 2 [1M]
## Fecha Pollo
## <mth> <dbl>
## 1 2001 Aug 65.6
## 2 2001 Sep 66.5
## 3 2001 Oct 65.7
## 4 2001 Nov 64.3
## 5 2001 Dec 63.2
## 6 2002 Jan 62.9
## 7 2002 Feb 62.9
## 8 2002 Mar 62.7
## 9 2002 Apr 62.5
## 10 2002 May 63.4
## # … with 170 more rows
## # ℹ Use `print(n = ...)` to see more rowsNote que se ajusta una tendencia que se ve que no es lineal.
tibble_chicken%>%timetk::plot_time_series(Fecha, Pollo,
.interactive = interactive,
.plotly_slider = TRUE)
###Usa Loess para hacer el ajuste de la tendencia, es decir usar smooth_vec() como versión simplificada de stats::loess()
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.75, degree = 2))## # A tibble: 180 × 3
## Fecha Pollo Pollo_ajus
## <date> <dbl> <dbl>
## 1 2001-08-01 65.6 62.8
## 2 2001-09-01 66.5 63.0
## 3 2001-10-01 65.7 63.2
## 4 2001-11-01 64.3 63.4
## 5 2001-12-01 63.2 63.6
## 6 2002-01-01 62.9 63.8
## 7 2002-02-01 62.9 64.0
## 8 2002-03-01 62.7 64.3
## 9 2002-04-01 62.5 64.5
## 10 2002-05-01 63.4 64.7
## # … with 170 more rows
## # ℹ Use `print(n = ...)` to see more rowstibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.5, degree = 1))%>%
ggplot(aes(Fecha, Pollo)) +
geom_line() +
geom_line(aes(y = Pollo_ajus), color = "red")
Descomposición
Vamos a usar la descomposición por medio de filtros de promedios móviles
chicken_decompo=decompose(chicken)
plot(chicken_decompo)
chicken_decompo$trend## Jan Feb Mar Apr May Jun Jul
## 2001
## 2002 NA 63.91958 63.75667 63.54250 63.32875 63.15500 63.05458
## 2003 63.64792 63.96792 64.36708 64.81542 65.31667 65.89792 66.51458
## 2004 71.92417 72.96750 73.82292 74.50458 75.07875 75.53583 75.88917
## 2005 75.35458 74.87458 74.52875 74.33750 74.18583 74.00208 73.75333
## 2006 71.02208 70.63875 70.27000 69.88542 69.53458 69.30333 69.28708
## 2007 73.70708 74.62875 75.53333 76.40667 77.19292 77.87083 78.43000
## 2008 80.89500 81.48792 82.07125 82.68125 83.38750 84.19292 85.03333
## 2009 87.24208 87.18625 86.97083 86.62875 86.23333 85.83042 85.45208
## 2010 84.68500 84.69750 84.85958 85.14083 85.44125 85.71333 85.92833
## 2011 86.39167 86.38500 86.45042 86.59625 86.84792 87.19125 87.60042
## 2012 91.02625 91.62125 92.18625 92.74917 93.34292 93.97792 94.66958
## 2013 99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167 NA NA NA NA NA NA
## Aug Sep Oct Nov Dec
## 2001 NA NA NA NA NA
## 2002 63.03542 63.09125 63.18125 63.27625 63.41958
## 2003 67.17167 67.90875 68.76083 69.72792 70.79583
## 2004 76.14000 76.26292 76.26458 76.13750 75.82708
## 2005 73.41375 72.99042 72.48750 71.95000 71.45333
## 2006 69.53958 70.06750 70.84500 71.77125 72.74708
## 2007 78.90083 79.32750 79.69417 80.02250 80.39333
## 2008 85.76458 86.26667 86.59333 86.87833 87.12667
## 2009 85.13500 84.92125 84.84500 84.81958 84.75667
## 2010 86.08375 86.24417 86.37750 86.42792 86.42208
## 2011 88.07750 88.61125 89.17625 89.77500 90.40333
## 2012 95.41000 96.14625 96.89542 97.70167 98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016Simule unos datos III, crear un objeto con periodicidad s=12, y extraer la descomposición, qué puede observar?
y=ts(rnorm(1000,0,1),start=c(2000,01),frequency = 4)
plot(decompose(y))
Descomposición STL
STL son las iniciales de “Seasonal and Trend decomposition using Loess”,el cual fue desarrollado por R. B. Cleveland et al. (1990).
library(feasts)
library(fable)
tsibble_chicken<-as_tsibble(chicken)
str(tsibble_chicken)## tbl_ts [180 × 2] (S3: tbl_ts/tbl_df/tbl/data.frame)
## $ index: mth [1:180] 2001 Aug, 2001 Sep, 2001 Oct, 2001 Nov, 2001 Dec, 2002 Jan...
## $ value: num [1:180] 65.6 66.5 65.7 64.3 63.2 ...
## - attr(*, "key")= tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
## ..$ .rows: list<int> [1:1]
## .. ..$ : int [1:180] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..@ ptype: int(0)
## - attr(*, "index")= chr "index"
## ..- attr(*, "ordered")= logi TRUE
## - attr(*, "index2")= chr "index"
## - attr(*, "interval")= interval [1:1] 1M
## ..@ .regular: logi TRUEtsibble_chicken %>%
model(
STL(value ~ trend() +
season(window = "periodic"),
robust = TRUE)) %>%
components() %>%
autoplot()
Note que en ambos casos se obliga a extraer un componente estacional, sin embargo puede que está componente en verdad no exista, por eso se debe verificar que en efecto hay.
set.seed(154)
w = rnorm(200); x = cumsum(w)
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)
La diferencia ordinaria de orden 1 es , \[\nabla^1 Y_t=(1-B)^1 Y_t=Y_t-Y_{t-1}\]
dx=diff(x)
plot.ts(dx, ylim=c(-5,5), main="Serie Diferenciada", ylab='')
par(mar = c(2,2,2,2))
fit = lm(chicken~time(chicken), na.action=NULL) # Regresión sobre el tiempo
par(mfrow=c(2,1))
plot(resid(fit), type="l", main="sin tendencia")
plot(diff(chicken), type="l", main="Primera Diferencia") 
par(mar = c(3,2,3,2))
par(mfrow=c(3,1)) # plot ACFs
acf(chicken, 48, main="ACF Pollo")
acf(resid(fit), 48, main="ACF Sin tendencia")
acf(diff(chicken), 48, main="ACF Primera Diferencia")
Transformación de Box-Cox para Estabilizar la Varianza Marginal
En ocasiones la serie presenta varianza marginal no constante a lo largo del tiempo, lo cual hace necesario tener en cuenta tal característica. En este caso, se siguiere hacer una transformación de potencia para estabilizar la varianza. Esta familia de transformaciones se llaman transformaciones Box-Cox.
\[ f_{\lambda}(u_{t})= \begin{cases} \lambda^{-1}(u^{\lambda}_{t}-1), & \text{si $u_{t} \geq 0$, para $\lambda>0$,}\\ \ln(u_{t}), &\text{ si $u_{t}>0$, para $\lambda=0$}. \end{cases} \]
data("AirPassengers")
plot(AirPassengers)
#####Transformación Box-Cox
library(FitAR)## Loading required package: lattice## Loading required package: leaps## Loading required package: ltsa## Loading required package: bestglmlibrary(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'forecast'## The following object is masked from 'package:FitAR':
##
## BoxCox## The following object is masked from 'package:astsa':
##
## gasforecast::BoxCox.lambda(AirPassengers, method = "guerrero", lower = -1, upper = 3) ###Me entrega el valor de lambda ## [1] -0.294731##method="loglik"
FitAR::BoxCox(AirPassengers)###Me entrega una gráfica
plot(forecast::BoxCox(AirPassengers,lambda=-0.294731))###
lAirPass=log(AirPassengers)
x11()
par(mar = c(1,1,1,1))
par(mfrow=c(2,1))
plot(AirPassengers,main="Serie de Pasajeros sin Transformar")
plot(lAirPass,main="Series con Transformación BoxCox")
##Box-Cox con timetk
timetk::box_cox_vec(AirPassengers,lambda = 'auto',silent = F)## box_cox_vec(): Using value for lambda: -0.294715585559316## Jan Feb Mar Apr May Jun Jul Aug
## 1949 2.548484 2.561375 2.588408 2.582937 2.567506 2.593720 2.615089 2.615089
## 1950 2.555038 2.577299 2.603899 2.593720 2.575381 2.616631 2.646225 2.646225
## 1951 2.610379 2.618160 2.656279 2.636912 2.648795 2.656279 2.680103 2.680103
## 1952 2.647515 2.658701 2.673640 2.659900 2.662270 2.699009 2.709884 2.720049
## 1953 2.676904 2.676904 2.715050 2.714201 2.709007 2.720866 2.737089 2.742835
## 1954 2.685298 2.668053 2.714201 2.707236 2.713347 2.737089 2.762580 2.756933
## 1955 2.720049 2.712489 2.739270 2.740706 2.741419 2.770363 2.796341 2.787869
## 1956 2.751056 2.746317 2.771524 2.769193 2.772100 2.801088 2.818144 2.814820
## 1957 2.770363 2.761963 2.792420 2.788382 2.791921 2.821786 2.837892 2.838594
## 1958 2.784223 2.772100 2.795371 2.788382 2.795857 2.826872 2.846724 2.851232
## 1959 2.794394 2.785275 2.815241 2.810978 2.820985 2.840332 2.864126 2.867216
## 1960 2.819775 2.808794 2.820583 2.836477 2.840332 2.860370 2.883509 2.879580
## Sep Oct Nov Dec
## 1949 2.595457 2.563441 2.529834 2.561375
## 1950 2.629937 2.590196 2.552878 2.602242
## 1951 2.663443 2.635540 2.611963 2.640966
## 1952 2.690331 2.671428 2.648795 2.674735
## 1953 2.715895 2.692301 2.658701 2.682201
## 1954 2.733382 2.709007 2.684272 2.709007
## 1955 2.768604 2.744238 2.715895 2.747003
## 1956 2.791921 2.765020 2.742129 2.765020
## 1957 2.814399 2.787869 2.764414 2.782096
## 1958 2.814399 2.793903 2.767420 2.782631
## 1959 2.837187 2.815659 2.795371 2.814820
## 1960 2.852177 2.836477 2.808352 2.825716Note que ahora usamos la misma función para verificar si en verdad la varianza fue estabilizada.
FitAR::BoxCox(lAirPass)
forecast::BoxCox.lambda(lAirPass, method = "guerrero", lower = -1, upper = 2)## [1] -0.419081library(VGAM)## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:tidyr':
##
## filllibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:VGAM':
##
## logit## The following object is masked from 'package:FitAR':
##
## Boot## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someVGAM::yeo.johnson(AirPassengers, lambda = 0)## Jan Feb Mar Apr May Jun Jul Aug
## 1949 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## 1950 4.753590 4.844187 4.955827 4.912655 4.836282 5.010635 5.141664 5.141664
## 1951 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## 1952 5.147494 5.198497 5.267858 5.204007 5.214936 5.389072 5.442418 5.493061
## 1953 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## 1954 5.323010 5.241747 5.463832 5.429346 5.459586 5.579730 5.713733 5.683580
## 1955 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## 1956 5.652489 5.627621 5.762051 5.749393 5.765191 5.926926 6.025866 6.006353
## 1957 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## 1958 5.831882 5.765191 5.894403 5.855072 5.897154 6.077642 6.198479 6.226537
## 1959 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## 1960 6.035481 5.971262 6.040255 6.135565 6.159095 6.284134 6.434547 6.408529
## Sep Oct Nov Dec
## 1949 4.919981 4.787492 4.653960 4.779123
## 1950 5.068904 4.897840 4.744932 4.948760
## 1951 5.220356 5.093750 4.990433 5.117994
## 1952 5.347108 5.257495 5.153292 5.273000
## 1953 5.472271 5.356586 5.198497 5.308268
## 1954 5.560682 5.438079 5.318120 5.438079
## 1955 5.746203 5.616771 5.472271 5.631212
## 1956 5.874931 5.726848 5.605802 5.726848
## 1957 6.003887 5.852202 5.723585 5.820083
## 1958 6.003887 5.886104 5.739793 5.823046
## 1959 6.139885 6.011267 5.894403 6.006353
## 1960 6.232448 6.135565 5.968708 6.070738car::yjPower(AirPassengers,lambda=0)## [1] 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## [9] 4.919981 4.787492 4.653960 4.779123 4.753590 4.844187 4.955827 4.912655
## [17] 4.836282 5.010635 5.141664 5.141664 5.068904 4.897840 4.744932 4.948760
## [25] 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## [33] 5.220356 5.093750 4.990433 5.117994 5.147494 5.198497 5.267858 5.204007
## [41] 5.214936 5.389072 5.442418 5.493061 5.347108 5.257495 5.153292 5.273000
## [49] 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## [57] 5.472271 5.356586 5.198497 5.308268 5.323010 5.241747 5.463832 5.429346
## [65] 5.459586 5.579730 5.713733 5.683580 5.560682 5.438079 5.318120 5.438079
## [73] 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## [81] 5.746203 5.616771 5.472271 5.631212 5.652489 5.627621 5.762051 5.749393
## [89] 5.765191 5.926926 6.025866 6.006353 5.874931 5.726848 5.605802 5.726848
## [97] 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## [105] 6.003887 5.852202 5.723585 5.820083 5.831882 5.765191 5.894403 5.855072
## [113] 5.897154 6.077642 6.198479 6.226537 6.003887 5.886104 5.739793 5.823046
## [121] 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## [129] 6.139885 6.011267 5.894403 6.006353 6.035481 5.971262 6.040255 6.135565
## [137] 6.159095 6.284134 6.434547 6.408529 6.232448 6.135565 5.968708 6.070738Ejercicio: Use la serie varve del paquete astsa para chequear si es necesario hacer transformación Box-Cox.
Vamos a correr lo mismo pero en Python para la transformación de BoxCox
import sys
print(sys.path)## ['', '/opt/anaconda3/envs/Python38andR/bin', '/opt/anaconda3/envs/Python38andR/lib/python38.zip', '/opt/anaconda3/envs/Python38andR/lib/python3.8', '/opt/anaconda3/envs/Python38andR/lib/python3.8/lib-dynload', '/Users/sergiocalderon/.local/lib/python3.8/site-packages', '/opt/anaconda3/envs/Python38andR/lib/python3.8/site-packages', '/Library/Frameworks/R.framework/Versions/4.1/Resources/library/reticulate/python']import pandas as pd
import matplotlib.pylab as plt
data = pd.read_csv('AirPassengers.csv')
print(data)## Month NPassengers
## 0 1949-01 112
## 1 1949-02 118
## 2 1949-03 132
## 3 1949-04 129
## 4 1949-05 121
## .. ... ...
## 139 1960-08 606
## 140 1960-09 508
## 141 1960-10 461
## 142 1960-11 390
## 143 1960-12 432
##
## [144 rows x 2 columns]print('\n Data Types:')##
## Data Types:print(data.dtypes)## Month object
## NPassengers int64
## dtype: objectcon=data['Month']
data['Month']=pd.to_datetime(data['Month'])
data## Month NPassengers
## 0 1949-01-01 112
## 1 1949-02-01 118
## 2 1949-03-01 132
## 3 1949-04-01 129
## 4 1949-05-01 121
## .. ... ...
## 139 1960-08-01 606
## 140 1960-09-01 508
## 141 1960-10-01 461
## 142 1960-11-01 390
## 143 1960-12-01 432
##
## [144 rows x 2 columns]pasajeros=data.set_index('Month')
#check datatype of index
#convert to time series:
ts = pasajeros['NPassengers']
ts.head(10)
####Graficar la Serie#####
#plt.plot(ts)
#plt.title('AirPassengers')
#plt.show()
## Month
## 1949-01-01 112
## 1949-02-01 118
## 1949-03-01 132
## 1949-04-01 129
## 1949-05-01 121
## 1949-06-01 135
## 1949-07-01 148
## 1949-08-01 148
## 1949-09-01 136
## 1949-10-01 119
## Name: NPassengers, dtype: int64print(r.AirPassengers)## [112.0, 118.0, 132.0, 129.0, 121.0, 135.0, 148.0, 148.0, 136.0, 119.0, 104.0, 118.0, 115.0, 126.0, 141.0, 135.0, 125.0, 149.0, 170.0, 170.0, 158.0, 133.0, 114.0, 140.0, 145.0, 150.0, 178.0, 163.0, 172.0, 178.0, 199.0, 199.0, 184.0, 162.0, 146.0, 166.0, 171.0, 180.0, 193.0, 181.0, 183.0, 218.0, 230.0, 242.0, 209.0, 191.0, 172.0, 194.0, 196.0, 196.0, 236.0, 235.0, 229.0, 243.0, 264.0, 272.0, 237.0, 211.0, 180.0, 201.0, 204.0, 188.0, 235.0, 227.0, 234.0, 264.0, 302.0, 293.0, 259.0, 229.0, 203.0, 229.0, 242.0, 233.0, 267.0, 269.0, 270.0, 315.0, 364.0, 347.0, 312.0, 274.0, 237.0, 278.0, 284.0, 277.0, 317.0, 313.0, 318.0, 374.0, 413.0, 405.0, 355.0, 306.0, 271.0, 306.0, 315.0, 301.0, 356.0, 348.0, 355.0, 422.0, 465.0, 467.0, 404.0, 347.0, 305.0, 336.0, 340.0, 318.0, 362.0, 348.0, 363.0, 435.0, 491.0, 505.0, 404.0, 359.0, 310.0, 337.0, 360.0, 342.0, 406.0, 396.0, 420.0, 472.0, 548.0, 559.0, 463.0, 407.0, 362.0, 405.0, 417.0, 391.0, 419.0, 461.0, 472.0, 535.0, 622.0, 606.0, 508.0, 461.0, 390.0, 432.0]print(r.lAirPass)## [4.718498871295094, 4.770684624465665, 4.882801922586371, 4.859812404361672, 4.795790545596741, 4.90527477843843, 4.997212273764115, 4.997212273764115, 4.912654885736052, 4.77912349311153, 4.6443908991413725, 4.770684624465665, 4.74493212836325, 4.836281906951478, 4.948759890378168, 4.90527477843843, 4.8283137373023015, 5.003946305945459, 5.135798437050262, 5.135798437050262, 5.062595033026967, 4.890349128221754, 4.736198448394496, 4.941642422609304, 4.976733742420574, 5.0106352940962555, 5.181783550292085, 5.093750200806762, 5.147494476813453, 5.181783550292085, 5.293304824724492, 5.293304824724492, 5.214935757608986, 5.087596335232384, 4.983606621708336, 5.111987788356544, 5.14166355650266, 5.19295685089021, 5.262690188904886, 5.198497031265826, 5.209486152841421, 5.384495062789089, 5.438079308923196, 5.488937726156687, 5.342334251964811, 5.25227342804663, 5.147494476813453, 5.267858159063328, 5.278114659230517, 5.278114659230517, 5.4638318050256105, 5.459585514144159, 5.43372200355424, 5.493061443340548, 5.575949103146316, 5.605802066295998, 5.4680601411351315, 5.351858133476067, 5.19295685089021, 5.303304908059076, 5.318119993844216, 5.236441962829949, 5.459585514144159, 5.424950017481403, 5.455321115357702, 5.575949103146316, 5.71042701737487, 5.680172609017068, 5.556828061699537, 5.43372200355424, 5.313205979041787, 5.43372200355424, 5.488937726156687, 5.4510384535657, 5.58724865840025, 5.594711379601839, 5.598421958998375, 5.752572638825633, 5.8971538676367405, 5.849324779946859, 5.7430031878094825, 5.6131281063880705, 5.4680601411351315, 5.627621113690637, 5.648974238161206, 5.6240175061873385, 5.75890177387728, 5.746203190540153, 5.762051382780177, 5.924255797414532, 6.023447592961033, 6.003887067106539, 5.872117789475416, 5.723585101952381, 5.602118820879701, 5.723585101952381, 5.752572638825633, 5.707110264748875, 5.87493073085203, 5.8522024797744745, 5.872117789475416, 6.045005314036012, 6.142037405587356, 6.1463292576688975, 6.0014148779611505, 5.849324779946859, 5.720311776607412, 5.817111159963204, 5.8289456176102075, 5.762051382780177, 5.8916442118257715, 5.8522024797744745, 5.8944028342648505, 6.075346031088684, 6.19644412779452, 6.22455842927536, 6.0014148779611505, 5.883322388488279, 5.736572297479192, 5.820082930352362, 5.886104031450156, 5.834810737062605, 6.0063531596017325, 5.981414211254481, 6.040254711277414, 6.156978985585555, 6.306275286948016, 6.326149473155099, 6.137727054086234, 6.008813185442595, 5.8916442118257715, 6.003887067106539, 6.0330862217988015, 5.968707559985366, 6.037870919922137, 6.133398042996649, 6.156978985585555, 6.282266746896006, 6.432940092739179, 6.406879986069314, 6.230481447578482, 6.133398042996649, 5.966146739123692, 6.068425588244111]py$data## Month NPassengers
## 1 1948-12-31 19:00:00 112
## 2 1949-01-31 19:00:00 118
## 3 1949-02-28 19:00:00 132
## 4 1949-03-31 19:00:00 129
## 5 1949-04-30 19:00:00 121
## 6 1949-05-31 19:00:00 135
## 7 1949-06-30 19:00:00 148
## 8 1949-07-31 19:00:00 148
## 9 1949-08-31 19:00:00 136
## 10 1949-09-30 19:00:00 119
## 11 1949-10-31 19:00:00 104
## 12 1949-11-30 19:00:00 118
## 13 1949-12-31 19:00:00 115
## 14 1950-01-31 19:00:00 126
## 15 1950-02-28 19:00:00 141
## 16 1950-03-31 19:00:00 135
## 17 1950-04-30 19:00:00 125
## 18 1950-05-31 19:00:00 149
## 19 1950-06-30 19:00:00 170
## 20 1950-07-31 19:00:00 170
## 21 1950-08-31 19:00:00 158
## 22 1950-09-30 19:00:00 133
## 23 1950-10-31 19:00:00 114
## 24 1950-11-30 19:00:00 140
## 25 1950-12-31 19:00:00 145
## 26 1951-01-31 19:00:00 150
## 27 1951-02-28 19:00:00 178
## 28 1951-03-31 19:00:00 163
## 29 1951-04-30 19:00:00 172
## 30 1951-05-31 19:00:00 178
## 31 1951-06-30 19:00:00 199
## 32 1951-07-31 19:00:00 199
## 33 1951-08-31 19:00:00 184
## 34 1951-09-30 19:00:00 162
## 35 1951-10-31 19:00:00 146
## 36 1951-11-30 19:00:00 166
## 37 1951-12-31 19:00:00 171
## 38 1952-01-31 19:00:00 180
## 39 1952-02-29 19:00:00 193
## 40 1952-03-31 19:00:00 181
## 41 1952-04-30 19:00:00 183
## 42 1952-05-31 19:00:00 218
## 43 1952-06-30 19:00:00 230
## 44 1952-07-31 19:00:00 242
## 45 1952-08-31 19:00:00 209
## 46 1952-09-30 19:00:00 191
## 47 1952-10-31 19:00:00 172
## 48 1952-11-30 19:00:00 194
## 49 1952-12-31 19:00:00 196
## 50 1953-01-31 19:00:00 196
## 51 1953-02-28 19:00:00 236
## 52 1953-03-31 19:00:00 235
## 53 1953-04-30 19:00:00 229
## 54 1953-05-31 19:00:00 243
## 55 1953-06-30 19:00:00 264
## 56 1953-07-31 19:00:00 272
## 57 1953-08-31 19:00:00 237
## 58 1953-09-30 19:00:00 211
## 59 1953-10-31 19:00:00 180
## 60 1953-11-30 19:00:00 201
## 61 1953-12-31 19:00:00 204
## 62 1954-01-31 19:00:00 188
## 63 1954-02-28 19:00:00 235
## 64 1954-03-31 19:00:00 227
## 65 1954-04-30 19:00:00 234
## 66 1954-05-31 19:00:00 264
## 67 1954-06-30 19:00:00 302
## 68 1954-07-31 19:00:00 293
## 69 1954-08-31 19:00:00 259
## 70 1954-09-30 19:00:00 229
## 71 1954-10-31 19:00:00 203
## 72 1954-11-30 19:00:00 229
## 73 1954-12-31 19:00:00 242
## 74 1955-01-31 19:00:00 233
## 75 1955-02-28 19:00:00 267
## 76 1955-03-31 19:00:00 269
## 77 1955-04-30 19:00:00 270
## 78 1955-05-31 19:00:00 315
## 79 1955-06-30 19:00:00 364
## 80 1955-07-31 19:00:00 347
## 81 1955-08-31 19:00:00 312
## 82 1955-09-30 19:00:00 274
## 83 1955-10-31 19:00:00 237
## 84 1955-11-30 19:00:00 278
## 85 1955-12-31 19:00:00 284
## 86 1956-01-31 19:00:00 277
## 87 1956-02-29 19:00:00 317
## 88 1956-03-31 19:00:00 313
## 89 1956-04-30 19:00:00 318
## 90 1956-05-31 19:00:00 374
## 91 1956-06-30 19:00:00 413
## 92 1956-07-31 19:00:00 405
## 93 1956-08-31 19:00:00 355
## 94 1956-09-30 19:00:00 306
## 95 1956-10-31 19:00:00 271
## 96 1956-11-30 19:00:00 306
## 97 1956-12-31 19:00:00 315
## 98 1957-01-31 19:00:00 301
## 99 1957-02-28 19:00:00 356
## 100 1957-03-31 19:00:00 348
## 101 1957-04-30 19:00:00 355
## 102 1957-05-31 19:00:00 422
## 103 1957-06-30 19:00:00 465
## 104 1957-07-31 19:00:00 467
## 105 1957-08-31 19:00:00 404
## 106 1957-09-30 19:00:00 347
## 107 1957-10-31 19:00:00 305
## 108 1957-11-30 19:00:00 336
## 109 1957-12-31 19:00:00 340
## 110 1958-01-31 19:00:00 318
## 111 1958-02-28 19:00:00 362
## 112 1958-03-31 19:00:00 348
## 113 1958-04-30 19:00:00 363
## 114 1958-05-31 19:00:00 435
## 115 1958-06-30 19:00:00 491
## 116 1958-07-31 19:00:00 505
## 117 1958-08-31 19:00:00 404
## 118 1958-09-30 19:00:00 359
## 119 1958-10-31 19:00:00 310
## 120 1958-11-30 19:00:00 337
## 121 1958-12-31 19:00:00 360
## 122 1959-01-31 19:00:00 342
## 123 1959-02-28 19:00:00 406
## 124 1959-03-31 19:00:00 396
## 125 1959-04-30 19:00:00 420
## 126 1959-05-31 19:00:00 472
## 127 1959-06-30 19:00:00 548
## 128 1959-07-31 19:00:00 559
## 129 1959-08-31 19:00:00 463
## 130 1959-09-30 19:00:00 407
## 131 1959-10-31 19:00:00 362
## 132 1959-11-30 19:00:00 405
## 133 1959-12-31 19:00:00 417
## 134 1960-01-31 19:00:00 391
## 135 1960-02-29 19:00:00 419
## 136 1960-03-31 19:00:00 461
## 137 1960-04-30 19:00:00 472
## 138 1960-05-31 19:00:00 535
## 139 1960-06-30 19:00:00 622
## 140 1960-07-31 19:00:00 606
## 141 1960-08-31 19:00:00 508
## 142 1960-09-30 19:00:00 461
## 143 1960-10-31 19:00:00 390
## 144 1960-11-30 19:00:00 432library(ggplot2)
ggplot2::ggplot(data = py$data,aes(x=Month,y=NPassengers) )+geom_line()
Gráficas de Retardos
Vamos a hacer gráficos de dispersión para chequear que tipos de relaciones hay entre los retardos de la variable interés. Vamos a trabajar con algunas series, por ejemplo: * Indice ambiental mensual (soi)(Southern Oscillation Index), el cual mide los cambios en la presión del aire, relacionados con las temperaturas de la superficie del mar en el Océano Pacífico central. * la serie rec (reclutamiento asociada al soi), número de nuevos peces. * Consumo mensual de gas natural en EE. UU.(USgas) medido en Billones de pies cúbicos. Esto permite chequear si hay posibles relaciones no-lineales.
library(astsa)
data("soi")
ts_info(soi)## The soi series is a ts object with 1 variable and 453 observations
## Frequency: 12
## Start time: 1950 1
## End time: 1987 9?soi
par(mar = c(2,2,2,2))
plot(soi, main="Indice soi")
par(mar = c(3,2,3,2))
astsa::lag1.plot(soi, 12) ###El 12 indica cuantos retardos y_t-k contra y_t 
###Hacer la gráfica con x11()En el gráfico de dispersión podemos ver que se muestra un ajuste no paramétrico, de la posible relación entre las variables al igual que una estimación de la autocorrelación entre \(s_t\) y \(s_{t-h}\). Vemos que varias de las relaciones exploradas parecen ser lineales. Uno puede observar que existen relaciones lineales positivas en los rezagos \(h = 1, 2, 11, 12\), mientras que negativas en los rezagos \(h=6,7\), las demás parecen ser no significativas o no lineales.
#pdf('/Users/sergiocalderon/Documents/Documentos - iMac de Sergio/Documentos iMac Sergio/Notas de Clase/Notas de clase/Notas de Clase Series de Tiempo Univariadas/Graficas/DispersionSoiRec.pdf',paper="USr")
par(mar = c(3,2,3,2))
lag2.plot(soi, rec, 8) #El 2 de lag2.plot es porque intervenienen dos serie de tiempo.
#dev.off()
lag2.plot(soi, rec, 8,corr=F)
Note también que el gráfico de dispersión del índice de nuevos peces con
retardos del indice soi nos muestra también relaciones posiblemente
lineales, al igual que no lineales. Ahora, podemos crear el gráfico de
autocorrelación simple que nos permite estimar y graficar las
autocorrelaciones para diferentes rezagos.
Note que con la función ts_lags de TSstudio podemos hacer una gráfica similar.
ts_lags(soi,lags=1:12)La función de autocorrelación simple o gráfico ACF y PACF
Cuando el proceso es estacionario, o al menos no presenta tendencia, podemos usar el gráfico acf para explorar las posibles relaciones lineales a diferentes rezagos. En seguida mostramos la función de autocorrelación para el índice soi y la serie de nuevos peces.
par(mfrow=c(2,1))
par(mar = c(2.7,2,2.7,2))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")
#dev.off()ccf(soi, rec, 48, main="SOI vs Recruitment", ylab="CCF")
#Corr(rec_{t},soi_{t-h})
#Corr(soi_{t},rec_{t-h},)index_soi_rec=as.Date(as.yearmon(tk_index(soi)))
df_soi_rec=data.frame(Fecha_soi_rec=index_soi_rec,soi=as.matrix(soi),rec=as.matrix(rec))
tibble_soi_rec=tibble(df_soi_rec)
tsibble::duplicates(tibble_soi_rec, key = NULL, index=Fecha_soi_rec) ##Mirar si hay registros duplicados## # A tibble: 0 × 3
## # … with 3 variables: Fecha_soi_rec <date>, soi <dbl>, rec <dbl>
## # ℹ Use `colnames()` to see all variable namestibble_soi_rec%>%plot_acf_diagnostics(Fecha_soi_rec,soi,.ccf_vars = rec,.lags = 36)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.En ambas gráficas podemos ver que se muestran periodicidades en las correlaciones que corresponden a valores separados por 12 unidades. También podemos ver que las observaciones con 12 meses o un año de diferencia están fuertemente correlacionadas positivamente, al igual que las observaciones en múltiplos como 24, 36, 48,\(\cdots\) Las observaciones separadas por seis meses están correlacionadas negativamente, lo que muestra que las excursiones positivas tienden a asociarse con las excursiones negativas a los seis meses ver Shumway2017 capítulo 1.
###AMI Del los libros H. Kantz and T. Schreiber: Nonlinear Time series Analysis (Cambridge university press) H. Abarbanel: Analysis of observed chaotic data (Springer, 1996) y NONLINEAR TIME SERIES ANALYSIS HOLGER KANTZ AND THOMAS SCHREIBER. Cambrige University Press 2003.
Ahora utilizaremos los paquetes nonlinearTseries y tseriesChaos para computar el average mutual information(AMI) o La información mutua promedio (AMI, la cual mide cuánto nos dice una variable aleatoria sobre otra, el cual se define como:
\[I(X;Y)=\sum_{i}\sum_{j}p(x_i,y_j)\log_2(\frac{p(x_i,y_j)}{p(x_i)p(y_j)}).\] En el contexto del análisis de series de tiempo, AMI ayuda a cuantificar la cantidad de conocimiento obtenido sobre el valor de \(X_{t+d}\) al observar \(X_t\). Equivalentemente, el AMI es una medida de qué tanto el conocimiento de \(X\) reduce la incertidumbre acerca de \(Y\). Esto implica que \(I(X,Y)=0\) si y sólo si \(X\) y \(Y\) son variables aletorias independientes. I(X; Y ) describe la información que la medición \(X_t\) en el tiempo \(t\) aporta a la medición \(X_{t+d}\) en el tiempo \(t + d\). Si se elige d como el valor alrededor del primer mínimo del AMI, entonces \(Y{t}\) e \(y_{t+d}\) son parcialmente pero no totalmente independientes.
Vamos a simular una serie de la forma \[x_t=\frac{x_t-1}{x_{t-3}^2+1}+\epsilon_t\] y a trabajar con la serie de linces Candienses del paquete astsa lynx.
library(nonlinearTseries)##
## Attaching package: 'nonlinearTseries'## The following object is masked from 'package:fabletools':
##
## estimate## The following object is masked from 'package:grDevices':
##
## contourLineslibrary(tseriesChaos)
et=rnorm(1100,0,1)
xt=rep(0,1100)
for(t in 13:1100)
{
xt[t]=(xt[t-1]-1)/(xt[t-12]^2+1)+et[t]
}
xtsimul=as.ts(xt[101:1100])
length(xtsimul)## [1] 1000plot(xtsimul)
acf(xtsimul)
par(mar = c(3,2,3,2))
astsa::lag1.plot(xtsimul, 6) 
nonlinearTseries::mutualInformation(xtsimul,lag.max = 100,n.partitions = 50,units = "Bits",do.plot = TRUE) #c("Nats", "Bits", "Bans")
## $time.lag
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [73] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
## [91] 90 91 92 93 94 95 96 97 98 99 100
##
## $mutual.information
## [1] 4.7211355 0.9086685 0.7868768 0.7146046 0.7234251 0.7393485 0.7257831
## [8] 0.7204692 0.7530212 0.7388565 0.6987539 0.7426234 0.8363724 0.7918314
## [15] 0.7495687 0.7384157 0.7227720 0.7306765 0.7185728 0.7490403 0.7888126
## [22] 0.7335948 0.7676440 0.7241823 0.7728012 0.7690717 0.7278799 0.7436484
## [29] 0.7654490 0.7561503 0.7324961 0.7820421 0.7413698 0.7871444 0.7796870
## [36] 0.7356812 0.7411850 0.7501907 0.7512091 0.7691925 0.7466313 0.7425748
## [43] 0.7152309 0.7567790 0.7658219 0.7516127 0.7371519 0.7288299 0.7622234
## [50] 0.7844343 0.7584129 0.7388515 0.7564102 0.7690676 0.7334833 0.7605943
## [57] 0.7808753 0.7568963 0.7571249 0.7135538 0.7925615 0.7584517 0.7768902
## [64] 0.7298138 0.7462323 0.7553250 0.7728154 0.7850680 0.7464423 0.7543733
## [71] 0.7465533 0.7813147 0.7620542 0.7528362 0.7663679 0.8111562 0.7237969
## [78] 0.7386509 0.7523882 0.7398758 0.7844187 0.7326143 0.7380267 0.7455960
## [85] 0.7691433 0.7309565 0.7479580 0.7314371 0.7583339 0.8090117 0.7704502
## [92] 0.7519965 0.7847387 0.7292170 0.7291928 0.7169975 0.7426216 0.7759557
## [99] 0.7898280 0.7688683 0.7700155
##
## $units
## [1] "Bits"
##
## $n.partitions
## [1] 50
##
## attr(,"class")
## [1] "mutualInf"pacf(xtsimul)
tseriesChaos::mutual(xtsimul, partitions = 50, lag.max = 100, plot=TRUE)
plot(astsa::Lynx)
acf(astsa::Lynx)
par(mar = c(3,2,3,2))
astsa::lag1.plot(astsa::Lynx, 12) 
tseriesChaos::mutual(astsa::Lynx, partitions = 50, lag.max = 10, plot=TRUE)
Detección de cíclos y estacionalidades
Podemos usar la función TSstudio::ts_heatmap para crear un mapa de calor. Este es un gráfico tridimensional, en donde en el eje y están los meses, y en el eje x están los años. Note en este caso que los meses de Diciembre,Enero, Febrero y Marzo presentan los valores mas oscuros, es decir los valores mas grandes a lo largo de los años, en contraste con los meses de Mayo a Septiembre que presentan colores mas claros. Este es un típico comportamiento de la presencia de un cíclo estacional en la serie. El flujo de color es horizontal.
TSstudio::USgas## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2000 2510.5 2330.7 2050.6 1783.3 1632.9 1513.1 1525.6 1653.1 1475.0 1567.8
## 2001 2677.0 2309.5 2246.6 1807.2 1522.4 1444.4 1598.1 1669.2 1494.1 1649.1
## 2002 2487.6 2242.4 2258.4 1881.0 1611.5 1591.4 1748.4 1725.7 1542.2 1645.9
## 2003 2700.5 2500.3 2197.9 1743.5 1514.7 1368.4 1600.5 1651.6 1428.6 1553.2
## 2004 2675.8 2511.1 2100.9 1745.2 1573.0 1483.7 1584.9 1578.0 1482.2 1557.2
## 2005 2561.9 2243.0 2205.8 1724.9 1522.6 1534.1 1686.6 1695.1 1422.5 1428.2
## 2006 2165.3 2144.4 2126.4 1681.0 1526.3 1550.9 1758.7 1751.7 1462.1 1644.2
## 2007 2475.6 2567.0 2128.8 1810.1 1559.1 1555.2 1659.9 1896.1 1590.5 1627.8
## 2008 2734.0 2503.4 2278.2 1823.9 1576.4 1604.2 1708.6 1682.9 1460.9 1635.8
## 2009 2729.7 2332.5 2170.7 1741.3 1504.0 1527.8 1658.0 1736.5 1575.0 1666.5
## 2010 2809.8 2481.0 2142.9 1691.8 1617.3 1649.5 1825.8 1878.9 1637.5 1664.9
## 2011 2888.6 2452.4 2230.5 1825.0 1667.4 1657.3 1890.5 1891.8 1655.6 1744.5
## 2012 2756.2 2500.7 2127.8 1953.1 1873.8 1868.4 2069.8 2008.8 1807.2 1901.1
## 2013 2878.8 2567.2 2521.1 1967.5 1752.5 1742.9 1926.3 1927.4 1767.0 1866.8
## 2014 3204.1 2741.2 2557.9 1961.7 1810.2 1745.4 1881.0 1933.1 1809.3 1912.8
## 2015 3115.0 2925.2 2591.3 2007.9 1858.1 1899.9 2067.7 2052.7 1901.3 1987.3
## 2016 3091.7 2652.3 2356.3 2083.8 1965.8 2000.7 2186.6 2208.4 1947.8 1925.2
## 2017 2914.2 2340.6 2523.7 1932.5 1892.5 1910.9 2142.1 2094.3 1920.9 2032.0
## 2018 3335.0 2705.9 2792.6 2346.3 2050.9 2058.7 2344.6 2307.7 2151.5 2279.1
## 2019 3399.9 2999.2 2899.9 2201.1 2121.0 2115.2 2407.5 2437.2 2215.6 2472.3
## Nov Dec
## 2000 1908.5 2587.5
## 2001 1701.0 2120.2
## 2002 1913.6 2378.9
## 2003 1753.6 2263.7
## 2004 1782.8 2327.7
## 2005 1663.4 2326.4
## 2006 1765.4 2122.8
## 2007 1834.5 2399.2
## 2008 1868.9 2399.7
## 2009 1776.2 2491.9
## 2010 1973.3 2714.1
## 2011 2031.9 2541.9
## 2012 2167.8 2503.9
## 2013 2316.9 2920.8
## 2014 2357.5 2679.2
## 2015 2249.1 2588.2
## 2016 2159.4 2866.3
## 2017 2357.7 3084.5
## 2018 2709.9 2993.1
## 2019 ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)TSstudio::ts_heatmap(USgas,title = "Mapa de Calor - Consumo de Gas Natural en EEUU")Voy enfocarme en un periodo de la serie de desempleo.
unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
title = "US Monthly Unemployment Rate",
Ytitle = "Unemployment Rate (%)",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)ts_info(unemployment)## The unemployment series is a ts object with 1 variable and 360 observations
## Frequency: 12
## Start time: 1990 1
## End time: 2019 12- El primer ciclo ocurrió entre 1990 y 2000, que estuvo cerca de un ciclo de 11 años.
- El segundo ciclo comenzó en 2000 y finalizó en 2007, que fue un ciclo de 7 años.
- Un tercer ciclo, que comenzó en 2007 y continúa a mayo de 2019 aún no se ha completado, lo que significa que se ha prolongado durante más de 12 años.
ts_heatmap(unemployment,title = "Mapa de Calor - Tasa de Desempleo EEUU Subserie")#ts_heatmap(USUnRate,title = "Mapa de Calor - Tasa de Desempleo EEUU")En este ejemplo, el flujo de color de la Tasa de Desempleo es vertical, lo que indica el estado del ciclo. En este caso, las franjas verticales más claras representan el final de un ciclo y el comienzo del siguiente. Asimismo, las franjas verticales más oscuras representan los picos del ciclo.
Explorando mas herramientas para detección de Estacionalidad
Medidas descriptivas
USgas_df <- data.frame(year = floor(time(USgas)), month = cycle(USgas),USgas = as.numeric(USgas))
USgas_df$month <- factor(month.abb[USgas_df$month], levels = month.abb)
library(dplyr)
USgas_summary <- USgas_df %>%group_by(month) %>%summarise(mean= mean(USgas),sd = sd(USgas))
USgas_summary## # A tibble: 12 × 3
## month mean sd
## <fct> <dbl> <dbl>
## 1 Jan 2806. 309.
## 2 Feb 2502. 223.
## 3 Mar 2325. 241.
## 4 Apr 1886. 175.
## 5 May 1708. 194.
## 6 Jun 1691. 216.
## 7 Jul 1864. 261.
## 8 Aug 1889. 237.
## 9 Sep 1687. 242.
## 10 Oct 1788. 260.
## 11 Nov 2015. 284.
## 12 Dec 2543. 278. library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout plot_ly (data = USgas_summary, x = ~ month, y = ~ mean, type = "bar", name = "Mean") %>%
layout (title = "USgas - Monthly Average", yaxis =list(title = "Mean", range = c(1500, 2700))) monthplot(USgas)
####También lo podemos hacer usando el paquete feasts y el objeto tsibble
USgas_tsbl=as_tsibble(USgas)
require(feasts)
USgas_tsbl=as_tsibble(USgas)
USgas_tsbl%>%gg_subseries(value)###Puede usar el argumento period=12 y da el mismo resultado, lo que significa es que se pueden agrupar las observaciones que están cada 12.
Note que basados en las estadísticas descriptivas podemos ver que las
medias son distintas para algunos meses, incluso sin quitar la
tendencia. A una misma conclusión llegamos basados en los gráficos.
Estas son típicas características de que hay presente un ciclo
estacional en la serie.
Explorando múltiples estacionalidades
Pueden haber múltiples estacionalidades(por lo general ocurren en series de alta frecuencia: diaria, cada hora, cada media y asi sucesivamente.) o una única estacionalidad, o inclusive cíclos que no se ven con facilidad.
Vamos a trabajar la base datos relacionada con el sistema nacional de transmisión de electricidad del Reino Unido. Más específicamente con la demanda de energía de forma horaria.
library(UKgrid)
require(TSstudio)
require(timetk)
require(feasts)
require(tsibble)
require(plotly)
UKgrid_xts <- extract_grid(type = "xts",
columns = "ND",
aggregate = "hourly",
na.rm = TRUE)
#extract_grid solo funciona para el conjunto de datos UKgrid
ts_plot(UKgrid_xts,
title = "National Hourly Demand UK Grid",
Ytitle = "Megawatts",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)UKgrid_tstbl <- extract_grid(type = "tsibble",
columns = "ND",
aggregate = "hourly",
na.rm = TRUE)
UKgrid_tbl <-as_tibble(UKgrid_tstbl)Como indicamos anteriormente, la primera indicación de la posible existencia de múltiples patrones estacionales en la serie es una frecuencia alta, como por diaria, horaria y en minutos. En esos casos, hay más de una forma de establecer la frecuencia de la serie. Por ejemplo, si capturamos una serie de tiempo de frecuencia diaria, la frecuencia de la serie se puede configurar de la siguiente manera: * Diariamente (o 365), asumiendo que el ciclo más apropiado es un año completo. * Días de semana (o 7) siempre que la oscilación del día de la semana sea más dominante que la del ciclo de año completo.
Al utilizar estadísticas descriptivas con este tipo de series, tendrá sentido aplicar este método para cada frecuencia potencial de la serie (o al menos las principales) con el fin de examinar si existe una indicación del patrón estacional.
Por ejemplo, UKgrid es una serie de tiempo por horas, que la marca automáticamente como sospechosa de tener múltiples patrones estacionales. Potencialmente, como se mencionó anteriormente, la demanda horaria de electricidad podría tener tres patrones estacionales diferentes:
- Horaria: Este es probablemente el principal patrón estacional de la serie, ya que existe una relación directa entre la demanda de energía eléctrica y la hora del día (hay alta demanda durante el día y baja demanda durante la noche o al contrario).
- Día de la semana: La demanda de electricidad a lo largo del día se deriva, potencialmente, del día de la semana. Tendría sentido esperar un alto consumo durante los días laborables y una menor tasa de consumo durante el fin de semana.
- Mensual: como los patrones climáticos varían a lo largo del año, la cantidad de luz del día y otros factores estacionales podrían afectar la demanda de electricidad.
Usando el paquete lubridate crearemos las características.
library(xts)
UKgrid_df <- data.frame(time = zoo::index(UKgrid_xts), UKgrid=as.numeric(UKgrid_xts))
str(UKgrid_df)## 'data.frame': 127296 obs. of 2 variables:
## $ time : POSIXct, format: "2005-04-01 00:00:00" "2005-04-01 01:00:00" ...
## $ UKgrid: num 65080 68207 69172 66769 64660 ...Ahora crearemos características estacionales basados en los periodos que deseamos explorar, por ejemplo hora del día,o día de la semana, o mes del año.
library(lubridate)
UKgrid_df$hour <- hour(UKgrid_df$time)
UKgrid_df$weekday <- wday(UKgrid_df$time, label = TRUE, abbr = TRUE)
UKgrid_df$month <- factor(month.abb[month(UKgrid_df$time)], levels = month.abb)
head(UKgrid_df)## time UKgrid hour weekday month
## 1 2005-04-01 00:00:00 65080 0 Fri Apr
## 2 2005-04-01 01:00:00 68207 1 Fri Apr
## 3 2005-04-01 02:00:00 69172 2 Fri Apr
## 4 2005-04-01 03:00:00 66769 3 Fri Apr
## 5 2005-04-01 04:00:00 64660 4 Fri Apr
## 6 2005-04-01 05:00:00 65209 5 Fri AprVamos a empezar las exploraciones analizando el ciclo horario.
UKgrid_hourly <- UKgrid_df %>%
dplyr::group_by(hour) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE), sd = sd(UKgrid, na.rm
= TRUE))
str(UKgrid_hourly)## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
## $ hour: int [1:24] 0 1 2 3 4 5 6 7 8 9 ...
## $ mean: num [1:24] 55622 54343 52920 51546 50509 ...
## $ sd : num [1:24] 9259 9482 9485 9165 8557 ...UKgrid_hourly## # A tibble: 24 × 3
## hour mean sd
## <int> <dbl> <dbl>
## 1 0 55622. 9259.
## 2 1 54343. 9482.
## 3 2 52920. 9485.
## 4 3 51546. 9165.
## 5 4 50509. 8557.
## 6 5 51747. 8657.
## 7 6 59048. 10489.
## 8 7 68627. 13347.
## 9 8 73581. 13003.
## 10 9 76320. 12557.
## # … with 14 more rows
## # ℹ Use `print(n = ...)` to see more rowsVamos ahora a hacer la gráfica de la media y la desviación estándar con base en las horas. Note que esas gráficas en escalas diferentes, así que hay que usar un gráfico especial.
require(plotly)
plot_ly(UKgrid_hourly) %>%
add_lines(x = ~ hour, y = ~ mean, name = "Media") %>%
add_lines(x = ~ hour, y = ~ sd, name = "Desviación Estándar", yaxis =
"y2",
line = list(color = "red", dash = "dash", width = 3)) %>%
layout(
title = "La demanda nacional de electricidad - Promedio horario vs. Desviación Estándar",
yaxis = list(title = "Media"),
yaxis2 = list(overlaying = "y",
side = "right",
title = "Desviación Estándar"
),
xaxis = list(title="Hora del Día"),
legend = list(x = 0.05, y = 0.9),
margin = list(l = 50, r = 50)
)Qué podemos destacar del gráfico y de las estadísticas descriptivas?
Hay poca demanda durante la noche (entre la medianoche y las 6 a.m.) y una alta demanda entre las horas de la mañana y la tarde.
Existe una fuerte correlación entre la demanda promedio y su desviación estándar.
La relativamente baja desviación estándar de la demanda promedio durante la noche podría indicar que existe un fuerte efecto subestacional durante esas horas además de la estacionalidad horaria. Esto debería tener sentido, ya que son horas de sueño normales y, por lo tanto, en promedio, la demanda es razonablemente la misma durante los días de semana.
Por otro lado, la alta desviación estándar a lo largo de las horas de alta demanda podría indicar que la demanda se distribuye de manera diferente en diferentes vistas de periodicidad (como día de la semana o mes del año).
Vamos a explorar el último punto, viendo la demanda a la madrugada(3 a.m) y empezando el día(9 a.m) con respecto al día de la semana.
UKgrid_weekday <- UKgrid_df %>%
dplyr::filter(hour == 3 | hour == 9) %>%
dplyr::group_by(hour, weekday) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
sd = sd(UKgrid, na.rm = TRUE))## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.UKgrid_weekday$hour <- factor(UKgrid_weekday$hour)
plot_ly(data = UKgrid_weekday, x = ~ weekday, y = ~ mean, type =
"bar",color = ~ hour) %>%
layout(title = "The Hourly Average Demand by Weekday",
yaxis = list(title = "Mean", range = c(30000, 75000)),
xaxis = list(title = "Weekday"))## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsEn el gráfico de barras anterior podemos ver que la demanda de electricidad a las 3 a.m. es relativamente estable durante todos los días de la semana, con una ligera diferencia entre el promedio durante los días laborables y los días del fin de semana (alrededor de un 2% diferente). Por otro lado, existe una diferencia entre la demanda del día laborable y el fin de semana a las 9 a.m. (es decir, la demanda del lunes es en promedio un 28% superior a la del domingo). Como era de esperar, esos resultados se alinearon con nuestras expectativas anteriores.
Ahora podemos aprovechar esos conocimientos para examinar si existe un patrón estacional mensual en la serie. Ahora seleccionaremos las mismas horas (3 a.m. y 9 a.m.); sin embargo, esta vez agruparemos estos datos por mes (en lugar de días laborables):
UKgrid_month <- UKgrid_df %>%
dplyr::filter(hour == 3 | hour == 9) %>%
dplyr::group_by(hour, month) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
sd = sd(UKgrid, na.rm = TRUE))## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.UKgrid_month$hour <- factor(UKgrid_month$hour)
plot_ly(data = UKgrid_month, x = ~ month, y = ~ mean, type = "bar",color =
~ hour) %>%
layout(title = "The Hourly Average Demand by Weekday",
yaxis = list(title = "Mean", range = c(30000, 75000)),
xaxis = list(title = "Month"))## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsPodemos ver en el gráfico de barras del resumen de agregación mensual que, en promedio, la demanda durante la noche (3 a.m.) y la mañana (9 a.m.) varía a lo largo de los meses del año. Además, hay un cambio significativo en la demanda durante la noche en comparación con la agregación entre semana. La variación de la serie de mes a mes indica la existencia de estacionalidad mensual en la serie.
Usando gráficos de densidades y Boxplot para explorar la estacionalidad.
Otro enfoque para analizar patrones estacionales en datos de series de tiempo es trazar la distribución de las unidades de frecuencia mediante el uso de histogramas o gráficos de densidad. Esto nos permitirá examinar si cada unidad de frecuencia tiene una distribución única que puede distinguirla del resto de unidades. Para esto usaremos el paquete ggplot2. Vamos empezar con la serie USgas.
require(timetk)
UKgrid_tbl## # A tibble: 127,296 × 2
## TIMESTAMP ND
## <dttm> <int>
## 1 2005-04-01 00:00:00 65080
## 2 2005-04-01 01:00:00 68207
## 3 2005-04-01 02:00:00 69172
## 4 2005-04-01 03:00:00 66769
## 5 2005-04-01 04:00:00 64660
## 6 2005-04-01 05:00:00 65209
## 7 2005-04-01 06:00:00 72376
## 8 2005-04-01 07:00:00 82679
## 9 2005-04-01 08:00:00 88706
## 10 2005-04-01 09:00:00 91108
## # … with 127,286 more rows
## # ℹ Use `print(n = ...)` to see more rowsUKgrid_tbl%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = ND,.feature_set = c("hour","wday.lbl","month.lbl"))## Warning: Removed 42 rows containing non-finite values (stat_boxplot).UKgrid_tstbl%>%gg_subseries(ND,period=24)
library(ggplot2)
ggplot(USgas_df, aes(x = USgas)) +
geom_density(aes(fill = month)) +
ggtitle("USgas - Kernel Density Estimates by Month") +
facet_grid(rows = vars(as.factor(month)))
Podemos ver algunos indicios de un patrón estacional en la serie, ya que las gráficas de densidad no se superponen entre sí (con la excepción de algunos meses consecutivos, como mayo y junio). Además, podemos ver que, durante algunos meses, la forma de las distribuciones es más plana con colas largas (principalmente durante los meses de invierno, noviembre, diciembre y enero). Sin embargo, no olvidemos el efecto de la tendencia o el crecimiento de un año a otro (como sabemos del capítulo anterior, la serie de gas de Estados Unidos tuvo una tendencia lineal desde el año 2010) ya que no la eliminamos de la serie. Repitamos este proceso; esta vez quitando la tendencia de la serie USgas antes de graficarla. Vamos a eliminarle la tendencia, posteriormente los explicaremos.
USgas_df$USgas_detrend <- USgas_df$USgas - decompose(USgas)$trend
ggplot(USgas_df, aes(x = USgas_detrend)) +
geom_density(aes(fill = month)) +
ggtitle("USgas - Estimación de la densidad vía Kernel por mes") +
facet_grid(rows = vars(as.factor(month)))## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.## Warning: Removed 12 rows containing non-finite values (stat_density).
Podemos ver que hay un comportamiento similar pero ahora las colas de
las densidades estimadas son mas cortas.Pueden hacer lo mismo pero en
vez de quitar la tendencia, hacemos una diferenciación.
En el caso de que la distribución de la mayoría de las unidades de frecuencia sea plana con una cola larga, podría ser una indicación de múltiples patrones estacionales en la serie. Regresemos a la serie UKgrid y tracemos las gráficas de densidad de 24 horas:
UKgrid_df$hour <- as.factor(UKgrid_df$hour)
ggplot(UKgrid_df, aes(x = UKgrid)) +
geom_density(aes(fill = hour)) +
ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
facet_grid(rows = vars(as.factor(hour)))## Warning: Removed 14 rows containing non-finite values (stat_density).
Como observamos anteriormente con las tablas de resúmenes estadísticos,
la distribución de la demanda neta de electricidad durante la noche es
relativamente estable (de ahí la distribución no plana con colas cortas
en contraposición a la distribución plana con cola larga durante el
día). Si ahora hacemos un subconjunto con una de las horas durante el
día y trazamos su distribución por el día de la semana, deberíamos
esperar una superposición durante la noche y poder distinguir entre la
distribución durante los días de la semana y el fin de semana, en
contraposición a solo el día de la semana.
Por ejemplo, la siguiente gráfica representa la distribución de la demanda a las 9 a.m. a lo largo de los días de la semana. Puede ver que la distribución durante los días de la semana se distingue de la del fin de semana:
UKgrid_df$weekday <- as.factor(UKgrid_df$weekday)
UKgrid_df %>% dplyr::filter(hour == 0) %>%
ggplot(aes(x = UKgrid)) +
geom_density(aes(fill = as.factor(weekday))) +
ggtitle("UKgrid - Kernel Density Estimates by Hour of the day") +
facet_grid(rows = vars(as.factor(weekday)))
Más herramientas del análisis de la estacionalidad
Vamos a usar el paquete forecast. Note que extraemos las subseries de los años. Lo cual demuestra un fuerte patrón estacional mensual.
library(forecast)
ggseasonplot(USgas,year.labels=TRUE,continuous=TRUE)
ggseasonplot(USgas, polar = TRUE)
Note que con el paquete TSstudio se puede hacer algo análogo.
ts_seasonal(USgas,type ="normal")ts_seasonal(USgas, type = "cycle")ts_seasonal(USgas, type = "box")ts_seasonal(USgas, type = "all")ts_seasonal con type=“cycle” añade un orden cronológico, en este caso por mes. Esto puede permitir la identificación de un patrón estacional sin tener que quitar la tendencia.
Cuando en el argumento type=“box”, elabora un gráfico de caja por unidad de frecuencia.
Finalmente cuando en el argumento type=“all”, elabora todos los gráficos anteriores.
El mapa de calor ya lo vimos anteriormente.
ts_heatmap(USgas, color = "Reds")Gráficos basados en cuantiles también son de utilidad.De forma predeterminada, la función devuelve un gráfico de cuantiles de las unidades de frecuencia de la serie, donde la línea media representa la mediana y las líneas inferior y superior representan los percentiles 25 y 75.
ts_quantile(UKgrid)## Warning in ts_quantile(UKgrid): The input object is a multiple time series
## object, by defualt will use only the first series as an inputEl argumento period permite examinar si los patrones estacionales de la serie están cambiando cuando se usa un subconjunto de tiempo diferente. Esto le permite examinar si la serie tiene patrones estacionales adicionales. Por ejemplo, podemos trazar el ciclo de 24 horas de la serie UKgrid por día de la semana estableciendo el argumento de período en días de la semana:
ts_quantile(UKgrid, period = "weekdays", n = 2)## Warning in ts_quantile(UKgrid, period = "weekdays", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## inputComo vimos anteriormente con las gráficas de densidad, la demanda de electricidad durante el día es relativamente mayor durante los días de semana en comparación con los fines de semana. De la misma manera, puede trazar el ciclo de 24 horas por mes:
ts_quantile(UKgrid, period = "monthly", n = 2)## Warning in ts_quantile(UKgrid, period = "monthly", n = 2): The input object is
## a multiple time series object, by defualt will use only the first series as an
## inputLa principal ventaja de las gráficas de cuantiles de múltiples períodos (por ejemplo, días de la semana, meses, etc.) sobre la gráfica de densidad que usamos anteriormente es que la primera representa todas las unidades de frecuencia (y, en el caso de UKgrid serie, con una frecuencia de 24 horas), mientras que el segundo representa una sola unidad de frecuencia (por ejemplo, la densidad de las 9 am durante los días de semana).
Usando Regresión para Descubrir un cíclo
Vamos a considerar la serie \[x_t =A\cos(2\pi\omega t+\varphi)+w_t,\] y una simulación de tamaño \(T=500\) para generar dos series, en donde las varianzas de los ruidos son diferentes en cada caso. Más específicamente \(\omega=\frac{1}{50}, A=2, \varphi=0.6\pi\) y \(\sigma_w=1,5.\)
cs = 2*cos(2*pi*1:500/50 + .6*pi); w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)
plot.ts(cs, main=expression(2*cos(2*pi*t/50+.6*pi)))
plot.ts(cs+w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,1)))
plot.ts(cs+5*w, main=expression(2*cos(2*pi*t/50+.6*pi) + N(0,25)))
Note que cuando la varianza del ruido es mas grande, es mas complicado
ver el cíclo.
Vamos a re-escribir el proceso anterior usando la relación trigonométrica \(\cos(\alpha \pm \beta) = \cos(\alpha) \cos(\beta) \mp \sin(\alpha) \sin(\beta),\) con lo cual:
\[A \cos(2\pi\omega t + \varphi) = \beta_1 \cos(2\pi\omega t) + \beta_2 \sin(2\pi\omega t)\] donde \(\beta_1=A\cos(\varphi)\) y \(\beta_2=-A\sin(\varphi)\), así que el proceso queda escrito de la siguiente manera:
\[x_t = \beta_1 \cos(2\pi t/50) + \beta_2 \sin(2\pi t/50) + w_t.\] asimiendo que \(\omega=1/50\) es conocido. La transformación nos permitió convertir la regresión no lineal(en los parámetros) en una lineal. Así podemos estimar \(\beta_1\) y \(\beta_2\) usando una regresión lineal. Note también que si tomamos \(\beta_1^2+\beta_2^2\), entonces tenemos \(A^2=\beta_1^2+\beta_2^2\), así que, un estimador de la amplitud al cuadrado es \(\hat{A}^2=\hat{\beta}_1^2+\hat{\beta}_2^2.\) Note que los valores verdaderos de \(\beta_1=\) y \(\beta_2=\)
beta1teor=2*cos(.6*pi)
beta2teor=-2*sin(.6*pi)
set.seed(90210) # fijamos una semilla para que nos dé el mismo resultado
x = 2*cos(2*pi*1:500/10 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)
z2 = sin(2*pi*1:500/50)
summary(fit <- lm(x~0+z1+z2)) # sin intercepto.##
## Call:
## lm(formula = x ~ 0 + z1 + z2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5372 -4.0634 -0.2505 3.3048 16.1621
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z1 -0.12622 0.34339 -0.368 0.713
## z2 -0.09283 0.34339 -0.270 0.787
##
## Residual standard error: 5.429 on 498 degrees of freedom
## Multiple R-squared: 0.0004179, Adjusted R-squared: -0.003597
## F-statistic: 0.1041 on 2 and 498 DF, p-value: 0.9012beta1teor## [1] -0.618034beta2teor## [1] -1.902113par(mar = c(2,2,2,2))
par(mfrow=c(2,1))
plot.ts(x)
plot.ts(x, col=8, ylab=expression(hat(x)))
lines(fitted(fit), col=2)
El periodograma
El periodograma es una herramienta que permite detectar el valor de la frecuencia \(\omega\) que en general es desconocida.
Exploración de Múltiples Cíclos y el Periodograma
Vamos a explorar una herremienta que sirve para describir funciones periódicas generales. Fourier mostró que toda función periódica puede representarse como una suma de funciones sinusoidales de distinta amplitud y frecuencia. La idea entonces es generalizar el análisis anterior para un ciclo a suma de funciones armónicas con distintas frecuencias.
Dada la serie de longitud \(T\), se denomina periodo básico o de Fourier, a las fracciones exactas del tamaño muestral, es decir \[s_j=\frac{T}{j},\ \ j=1,2,\cdots,T/2.\]
El valor máximo de periodo básico se obtiene para \(j=1\) y es \(T\), es decir, observamos la onda una sola vez. El mínimo se obtiene para \(j=T/2\) y es 2 ya que no se puede observar periodos que duren menos de 2 observaciones. Como antes, suele usarse las frecuencias en lugar de los periodos para el ajuste de los ciclos, así que las frecuencias básicas o de Fourier se define como \[f_j=\frac{j}{T}\ \ j=1,2,\cdots,T/2.\]
Así, \(1/T\leq f_j\leq 1/2\), y el valor máximo sería \(f_j=0.5\) y se conoce como la frecuencia Nyquist. Con esto, una serie de tiempo \(Z_t\) que es periódica se puede representar mediante
\[Z_t=\mu+\sum_{j=1}^{T/2}A_j\sin(\omega_j t)+ \sum_{j=1}^{T/2}B_j\cos(\omega_j t)+a_t.\]
Note que éste modelo tiene tantos parámetros como observaciones, con
lo cual hay que buscar un procedimiento para seleccionar las frecuencias
que deben ser incluidas para explicar la evolución de la serie. Para
esto se utilizará la herramienta conocida como el periodograma. Se puede
verificar que la contribución de una onda a la varianza es dada por
\(\frac{\hat{R}^2}{2}\), así que ondas
asociadas con amplitudes grandes serán importantes en la explicación de
la variabilidad de la serie, mientras que ondas asociadas con amplitudes
bajas serán poco importantes. De forma análoga a como se mostró en la
sección anterior, los estimadores de los coeficientes \(A_j\) y \(B_j\) están dados por \[\hat{A}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\sen(\omega_j
t)\]
\[\hat{B}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\cos(\omega_j
t)\] y \[\hat{R}_j=\hat{A}_j^2+\hat{B}_j^2,\] donde
\(\omega_j=2\pi f_j\) y \(f=\frac{j}{T}\), con lo cual \[TS_z^2=\frac{T}{2}\sum_{j=1}^{T/2}\hat{R}_j^2.\]
Se le da el nombre de periodograma a la representación de cada \(\frac{T\hat{R}_j^2}{2}\) en función de la frecuencia \(\omega_j\) o \(f_j\). Así \[\begin{equation} \label{periodograma} I(f_j)=\frac{T\hat{R}_j^2}{2},\ \ \text{con}\ 1/T\leq f_j\leq 1/2 \end{equation}\] y \[\bar{I}=\sum_{j=1}^{T/2}\hat{R}_j^2=S_z^2.\]
Nos enfocaremos en las frecuencias básicas únicamente, lo cual no es restrictivo si \(T\) es muy grande, puesto que el número de frecuencias básicas será muy grande y siempre existirá alguna frecuencia básica muy próxima a la que puede interesarnos. El periodograma puede verse como una herramienta para la detección de posibles ciclos(ocultos) deterministas en una serie temporal. Por ejemplo, en una serie mensual estacional de periodo \(s=12\), esperamos encontrar un valor alto de periodograma para \(f=1/12\), pero también para \(f=j/12\), es decir, \(1/6,1/4,1/3\) que son armónicos del periodo estacional. Por otro lado, la serie puede tener otros ciclos no necesariamente ligado al periodo estacional, siendo el periodograma una buena herramienta para detectar estos posibles componentes. \
Podemos considerar la amplitud calculada para la frecuencia \(f_j\) como un promedio de las amplitudes existentes en las frecuencias situadas en el intervalo \(f_j\pm \frac{1}{2}T\). Con esto, se puede obtener el periodograma suavizado construyendo rectángulos con centro \(f_j\), base igual a \(1/T\) y alturas \(I(f_j)=T\hat{R}_j/2\). Otra forma de suavizar el periodograma es \[I(f)=\sum_{f_i-q}^{fi+q}p_iI(f_i)\ \ 0\leq f \leq 0.5,\]
donde \(q\) representa la ventana usada y los \(p_i\) son los pesos simétricos y no modifica el valor de la varianza \(S^2_z\) que es área bajo la curva.
Vamos a retomar el ejemplo anterior donde simulamos una serie de tamaño \(T=500\) con frecuencia \(\omega=1/50\).
plot.ts(x)
spectrum(x,log='no')
abline(v=1/10, lty=2,col="red")
###Otra forma de computarlo es vía la transformada de Fourier
Px = Mod(fft(x))^2
Freq=0:499/500
plot(Freq, Px, type='l')
u = which.max(Px[1:250])
sprintf("El valor de la frecuencia donde se máximiza el periodograma es %s",Freq[u])## [1] "El valor de la frecuencia donde se máximiza el periodograma es 0.1"Note que en cualquiera de las dos formas de obtener el periodograma, hemos descubierto la frecuencia \(\omega=1/50\) como frecuencia príncipal.
En ocasiones es necesario suavizar el periodograma porque encontramos varios picos que en verdad no son valores grandes.
spectrum(x,log='no',span=5)
#span vector of odd integers giving the widths of modified Daniell smoothers to be used to smooth the periodogram.
spectrum(x,log='no',span=c(5,5))
spectrum(x,log='no',span=c(2,2))
Qué sucede si no consideramos un ciclo determinístico a través de
funciones sinusoidales sino que lo consideramos a través de variables
Dummy? Es decir por ejemplo una componente estacional de periodo 4,
modelado a través de variables dummy. \[Y_t=\delta_1\gamma_{1,t}+\delta_2\gamma_{2,t}+\delta_3\gamma_{3,t}+\delta_4\gamma_{4,t}+\epsilon_{t}\]
donde las variables dummy \(\gamma_{i,t},
i=1,2,3,4\) indican si se está en el trimestre \(i\). Puede en periodograma detactar este
tipo de estacionalidades obtenidas a través de dummy estacionales?
library(dplyr)
library(rstanarm)## Loading required package: Rcpp## This is rstanarm version 2.21.1## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())##
## Attaching package: 'rstanarm'## The following object is masked from 'package:car':
##
## logit## The following objects are masked from 'package:VGAM':
##
## cauchy, dirichlet, exponential, laplace, logitlibrary(rstan)## Loading required package: StanHeaders## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:tidyr':
##
## extractlibrary(ggplot2)
library(bayesplot)## This is bayesplot version 1.9.0## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settings=4
deltas=seq(1:s)
####Ejemplo Fácil Dummy+ruido
T=200*s
epsilon=rnorm(T,0,1)
l=diag(s)
X=matrix(rep(t(l),T/s),ncol=ncol(l),byrow=TRUE)
Y=X%*%deltas+epsilon
plot(as.ts(Y))
acf(as.ts(Y))
#simul=data.frame(resp=Y,diseno=X)
#salida_dummy_ruido=stan_glm(resp~. -1,data=simul)
PeriodgramaY=spectrum(Y,log='no')
ubicacionY=which.max(PeriodgramaY$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para la serie es: %s",PeriodgramaY$freq[ubicacionY])## [1] "El valor de la frecuencia donde se máximiza el periodograma para la serie es: 0.25"sprintf("El periodo correspondiente es aporximadamente: %s",1/PeriodgramaY$freq[ubicacionY])## [1] "El periodo correspondiente es aporximadamente: 4"Veamos ahora con una serie de tiempo real. Usaremos la serie SOi y Recruitment para ilustrar su uso. Recordemos que ambas series son mensuales, así que las frecuencias que serán dibujadas en el eje x vendrán en múltiplos de \(\frac{1}{12}\). Usaremos ahora la función mvspec del paquete astsa.
library(astsa)
data(soi)
soi.per = astsa::mvspec(soi, log="no")
ubicacionsoi=which.max(soi.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para SOI es: %s",soi.per$freq[ubicacionsoi])## [1] "El valor de la frecuencia donde se máximiza el periodograma para SOI es: 1"n_soi <- length(soi.per$spec)
valor_seg_per=sort(soi.per$spec,partial=n_soi-1)[n_soi-1]
ubica_segundo_soi=which(soi.per$spec==valor_seg_per)
sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: %s",soi.per$freq[ubica_segundo_soi])## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para SOI es: 0.2"abline(v=soi.per$freq[ubicacionsoi], lty=2,col="blue")
abline(v=soi.per$freq[ubica_segundo_soi], lty=2,col="blue")
Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\),
el cual corresponde al obvio ciclo anual (periodo de 12 meses). La
segunda corresponde a \(\omega=\frac{1}{5}
\cdot \frac{1}{12}=\frac{1}{60}\), el cual corresponde a un
posible ciclo de periodo \(1/(1/60)=
60\) meses, es decir 5 años, lo cual puede ser debido al fenómeno
de El Niño.
data("rec")
rec.per = astsa::mvspec(rec, log="no")
ubicacionrec=which.max(rec.per$spec)
sprintf("El valor de la frecuencia donde se máximiza el periodograma para REC es: %s",rec.per$freq[ubicacionrec])## [1] "El valor de la frecuencia donde se máximiza el periodograma para REC es: 1"n_rec <- length(rec.per$spec)
valor_seg_per_rec=sort(rec.per$spec,partial=n_rec-1)[n_rec-1]
ubica_segundo_rec=which(rec.per$spec==valor_seg_per_rec)
sprintf("El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: %s",rec.per$freq[ubica_segundo_rec])## [1] "El valor de la frecuencia donde se alcanza el segundo máximo para el periodograma para REC es: 0.25"abline(v=rec.per$freq[ubicacionrec], lty=2,col="blue")
abline(v=rec.per$freq[ubica_segundo_rec], lty=2,col="blue")
Note que la primera frecuencia mas alta se encuentra en \(\omega=1\cdot \frac{1}{12}=\frac{1}{12}\),
el cual corresponde al obvio ciclo anual (periodo de 12 meses), igual
que para el casi del indice soi. La segunda corresponde a \(\omega=\frac{1}{4} \cdot
\frac{1}{12}=\frac{1}{48}\), el cual corresponde a un posible
ciclo de periodo \(1/(1/48)= 48\)
meses, es decir 4 años. lo cual puede ser debido al fenómeno de El Niño.
Esta actividad de banda ancha sugiere que el posible ciclo de El Niño es
irregular, pero tiende a rondar los cuatro años en promedio.
Tarea : Hacer lo mismo pero ahora usar la serie de pasajeros diferenciada, al igual que la serie UKgrid, mas específicamente la serie UKgrid_xts y la serie taylor del paquete forecast.
library(sarima) ##
## Attaching package: 'sarima'## The following object is masked from 'package:rstanarm':
##
## se## The following object is masked from 'package:astsa':
##
## sarima## The following object is masked from 'package:stats':
##
## spectrumx <- ts(sarima::sim_sarima(n=144, model = list(siorder=1,
nseasons=12, sigma2 = 1)),frequency = 12 )#(1-B)^{12} X_t = e_t
plot(x)
acf(x)
monthplot(x)
ts_seasonal(x, type = "all")spectrum(x,log='no',span=5)## Estimated spectral density
## series: x
## method: Smoothed Periodogram
## nseasons: 12
## frequency range: (-6,6]
##
## sort method for the table: decreasing magnitudes
##
## freq spec % Total Cum. %
## [1,] 1.0000000 0.5488446 0.09121551 0.09121551
## [2,] 0.9166667 0.5259173 0.08740510 0.17862061
## [3,] 1.0833333 0.5008903 0.08324572 0.26186633
## [4,] 0.8333333 0.3152458 0.05239244 0.31425877
## [5,] 2.0833333 0.3137395 0.05214210 0.36640087
## [6,] 2.0000000 0.3113759 0.05174928 0.41815014
astsa::mvspec(x, log="no")
library(sarima)
x <- ts(sim_sarima(n=288,model=list(sar=0.3, nseasons=12, sigma2 = 1)), frequency = 12 )#SAR(1)
plot(x)
acf(x)
monthplot(x)
ts_seasonal(x, type = "all")spectrum(x,log='no',span=5)## Estimated spectral density
## series: x
## method: Smoothed Periodogram
## nseasons: 12
## frequency range: (-6,6]
##
## sort method for the table: decreasing magnitudes
##
## freq spec % Total Cum. %
## [1,] 2.000000 0.2843993 0.02370651 0.02370651
## [2,] 2.041667 0.2799563 0.02333616 0.04704267
## [3,] 4.166667 0.2643718 0.02203709 0.06907975
## [4,] 4.083333 0.2621709 0.02185363 0.09093338
## [5,] 4.125000 0.2549984 0.02125575 0.11218914
## [6,] 1.958333 0.2365108 0.01971469 0.13190383
astsa::mvspec(x, log="no")
raices=polyroot(c(1,rep(0,11),-0.3))
raices## [1] 0.5527684+0.9574230i -0.9574230+0.5527684i -0.5527684-0.9574230i
## [4] 0.9574230-0.5527684i 0.0000000+1.1055369i -1.1055369-0.0000000i
## [7] 0.0000000-1.1055369i 1.1055369+0.0000000i -0.5527684+0.9574230i
## [10] -0.9574230-0.5527684i 0.5527684-0.9574230i 0.9574230+0.5527684iconj_raices=Conj(raices)
raices[12]*conj_raices[12]## [1] 1.222212+0iSuavizamiento en el contexto de Series de Tiempo
Promedio Móvil
El promedio móvil es un método es útil para descubrir ciertos rasgos en una serie de tiempo, como tendencias a largo plazo y componentes estacionales. En particular, si \(x_t\) representa las observaciones, entonces una forma de predecir predecir o estimar la tendencia de la serie es:
\[m_t=\sum_{j=-k}^{k}a_jx_{t-j},\] donde si\(a_j=a_{-j}\ge0\) y \(\sum_{j=-k}^{k}a_j=1\) se conoce como el promedio móvil simétrico de los datos.
# Filtro Boxcar
wgts = c(.5, rep(1,11), .5)/12 ###Los pesos de los filtros
soif = stats::filter(soi, sides=2, filter=wgts)
plot(soi)
lines(soif, lwd=2, col=4)
###Forma del Filtro
nwgts = c(rep(0,20), wgts, rep(0,20))
plot(nwgts, type="l", ylim = c(-.02,.1), xaxt='n', yaxt='n', ann=FALSE)
## Suavizamiento Kernel El suavizamiento kernel es un suavizador de
promedio móvil que utiliza una función de ponderación, o kernel, para
promediar las observaciones.Veamos ahora como queda el promedio
móvil:
\[m_t=\sum_{i=1}^{n}w_i(t)x_i\] donde \[w_i(t)=K(\frac{t-i}{b})/\sum_{j=1}^{n}K(\frac{t-j}{b})\] son los pesos y \(K(\cdot)\) es una función kernel. Este estimador es llamado el Estimador de Nadaraya-Watson. Usaremos la función ksmooth
plot(soi)
lines(ksmooth(time(soi), soi, "normal", bandwidth=1), lwd=2, col=4)
#par(fig = c(.65, 1, .65, 1), new = TRUE) # the insert
gauss = function(x) { 1/sqrt(2*pi) * exp(-(x^2)/2) }
x = seq(from = -3, to = 3, by = 0.001)
plot(x, gauss(x), type ="l", ylim=c(-.02,.45), xaxt='n', yaxt='n', ann=FALSE)
## Lowess Otro enfoque para suavizar un gráfico de tiempo es la
regresión del vecino más cercano. La técnica se basa en la regresión de
k vecinos más cercanos, en la que uno usa solo los datos \(\{x_{t−k/2}, ..., x_t, ..., x_{t+k/2}\}\)
para predecir \(x_t\) mediante
regresión, y luego establece \(m_t =
\hat{x}_t\). Primero, una cierta proporción de vecinos más
cercanos a \(x_t\) se incluyen en un
esquema de ponderación; los valores más cercanos a \(x_t\) en el tiempo obtienen más peso.
Luego, se utiliza una regresión ponderada robusta para predecir \(x_t\) y obtener los valores suavizados
\(mt\). Cuanto mayor sea la fracción de
vecinos más cercanos incluidos, más suave será el ajuste. El R se usa la
función lowess.
plot(soi)
lines(lowess(soi, f=.05), lwd=2, col=4) # Ciclo El Niño
lines(lowess(soi), lty=2, lwd=2, col=2) # tendencia (con span por defecto)
Suavizamiento Splines
Una forma obvia de suavizar los datos sería ajustar una regresión polinomial en términos del tiempo. Por ejemplo, un polinomio cúbico tendría \(x_t = m_t + w_t\) donde \(m_t =\beta_0 + \beta_1t + \beta_2t^2 + \beta_3t^3\). Entonces podríamos ajustar \(m_t\) mediante mínimos cuadrados ordinarios.
Una extensión de la regresión polinomial es dividir primero el tiempo \(t = 1,. . . , n\), en k intervalos, \([t_0 = 1, t_1]\), \([t_1 + 1, t_2],\cdots,\) \([t_{k − 1} + 1, t_k = n]\); los valores \(t_0\), \(t_1\), …, \(t_k\) se llaman nodos. Luego, en cada intervalo, se ajusta una regresión polinomial, normalmente de orden 3, y esto se llama splines cúbicos. Un método relacionado es suavizar splines, que minimiza el compromiso entre el ajuste y el grado de suavidad dado por
\[\sum_{t=1}^{n}[x_t-m_t]^2+\lambda\int(m_t'')^2 dt,\] donde \(m_t\) es un spline cúbico con nodos en cada tiempo t y el grado de suavidad es controlado por \(\lambda>0.\) El parámetro de suavizado en R es controlado por el argumento spar de la función smooth.spline del paquete stats.
plot(soi)
lines(smooth.spline(time(soi), soi, spar=.5), lwd=2, col=4)
lines(smooth.spline(time(soi), soi, spar= 1), lty=2, lwd=2, col=2)
La diferencia entre uan y otra gráfica es es grado de suavizamiento.
<> Tarea: Hacer lo mismo para la serie Recruitment.
#```{r} #library(“readr”) #library(zoo)
#write_csv(data.frame(soi=as.matrix(soi),date=as.Date(as.yearmon(time(soi)))), file = “soi.csv”)
#write_csv(data.frame(rec=as.matrix(rec),date=as.Date(as.yearmon(time(rec)))), file = “rec.csv”) #```